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Using a combined analytical/molecular dynamics (MD) approach, we study the current fluctuation 
spectra and longitudinal and transverse collective mode dispersions of the classical two-dimensional 
(point) dipole system (2DDS) characterized by the 4>D{r) = ji^ jr'^ repulsive interaction potential; 
[1, is the electric dipole strength. The interest in such two-dimensional dipole systems (2DDS) is 
twofold. First, the quasi-long range interaction makes the system a unique classical many body 
system, with a remarkable collective mode behavior. Second, the system may be a good model for 
a closely spaced semiconductor electron-hole bilayer, a system that is in the forefront of current 
experimental interest. The longitudinal collective excitations, which are of primary interest for the 
liquid phase, are acoustic at long wavelengths. At higher wave numbers and for sufficiently high 
coupling strength, we observe the formation of a deep minimum in the dispersion curve preceded by a 
sharp maximum; this is identical to what has been observed in the dispersion of the zero temperature 
bosonic dipole system, which in turn emulates so called roton-maxon excitation spectrum of the 
superfluid ''He. The analysis we present gives an insight into the emergence of this apparently 
universal structure, governed by strong correlations. We study both the liquid and the crystalline 
solid state. We also observe the excitation of combination frequencies, resembling the roton-roton, 
roton-maxon, etc. structures in ''He. 

PACS numbers: 52.27.Gr,52.65.Yy,73.20.Mf 



I. INTRODUCTION 

This paper addresses the dynamical behavior of two- 
dimensional classical liquid and solid systems consisting 
of point electric dipoles in a plane, with the dipole mo- 
ments oriented perpendicular to the plane. This system 
is equivalent to an ensemble of point particles with a re- 
pulsive Xjr^ interaction potential. The interest in such 
two-dimensional dipole systems (2DDS) is twofold. First, 
the quasi-long range l/r^ interaction makes the collective 
behavior of the system somewhat similar to that of the 
one-component plasma with its l/r Coulomb interaction, 
while, at the same time, the Xjr^ interaction, with its 
rapid drop-off at large r and its hard singularity at r = 0, 
resembles the typical interaction potentials in classical 
liquids. This duality makes the 2DDS a unique classical 
many body system worthy of detailed exploration. The 
second aspect that makes the study of the 2DDS timely 
and important is that it can provide, as we will discuss 
below, a good model for a semiconductor electron-hole 
bilayer, a system that is in the forefront of current ex- 
perimental interest jl| , as well as for colloidal suspensions 
of super-paramagnetic particles [1] . 

The existence of bound electron-hole excitons in semi- 
conductors was predicted quite some time ago by Keldysh 
and co-workers [3, |3| and by Halperin and Rice [5|. 
Electron-hole bilayers (EHBs) are structures especially 



well suited to the formation of a stable dipole-like exci- 
tonic phase In such systems, the charges in the 

two layers with opposite polarities are physically sep- 
arated from each other, reducing their recombination 
rate and for sufficiently small layer separations forming 
a bound dipole-like excitonic structure. The formation 
of the excitonic/dipole phase in the EHB has been con- 
firmed by recent diffusion and path integral Monte Carlo 
(MC) simulations [ol-il l| . by classical MC simulations 
and by classical molecular dynamics (MD) simu- 
lations [l^l- Fixed- node diffusion MC simulations of the 
zero-temperature symmetric (me = mh, Ue = nu) EHB 
indicate three phases: the excitonic liquid, the spin- 
unpolarized Coulomb liquid, and the triangular Wigner 
crystal. Classical MC simulations of the bipolar bilayer 
\V2\ indicate the existence of four phases in the strong 
coupling regime: the excitonic dipole liquid and solid and 
the Coulomb liquid and solid phases. The necessary ex- 
istence of these four phases was also pointed out in [T^] ■ 
Since the excitons are bosons, at low temperatures they 
may also form a Bose-Einstein condensate fisl- Tsj , be- 
come superfluid pj'l , or possibly a supersolid [14 1 . 

Based on the phase diagrams the closely spaced 

EHB in its dipole-like excitonic phase can be, in a good 
approximation, modeled as a 2D monolayer of interact- 
ing point electric dipoles. The model 2D dipole system 
(2DDS) can be described as a collection of TV spinless 
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point dipoles, each of mass m = me + rrih, occupying the 
large but bounded area A; n = N/A is the average den- 
sity. The dipoles are free to move in the x — y plane with 
dipolar moment oriented in the z-direction; the repulsive 
interaction potential is accordingly given by (fio — /i^/r'^, 
where /i is the electric dipole strength. 

The coupling strength of the symmetric (rie = Uh = n, 
nie = mfi —771/2) EHB at arbitrary degeneracy is char- 
acterized by the parameter F — e^/(a(i?kin)), where e 
is the electrical charge and a = ^JYJ(^Tn) is the 2D 
Wigner-Seitz radius. In the high-temperature classical 
domain, this becomes the customary coupling parameter 
r = /3e^/a {P = 1/k-BT), while at zero-temperature, it 
becomes — a/ae, (ob — f^jm^^he^ is the Bohr ra- 
dius). By the same token, the high-temperature classi- 
cal coupling parameter for the 2DDS is characterized by 
Vu = whereas at zero-temperature, the appro- 

priate measure of the 2DDS coupling strength is given by 
= "To/a, where rp = m/i^/?i^ is the dipole equivalent 
of the Bohr radius; note the correspondence Tjj <^ r^j. 

In the case of the symmetric EHB in the zero- 
temperature quantum domain, the Coulombic and dipole 
coupling parameters are related to each other by r/j = 
2rs{d/a)'^, where d is the spacing between layers. High 
coupling (r£) ^ 1) for point dipoles corresponds to the 
low-density regime in the closely spaced {d/a < 1) EHB, 
as dictated by the ordering (2rs)~^ ^ {d/aY < 1, or 
equivalently, a > d ^ 03/2. In the high-temperature 
classical domain, the Coulombic and dipole coupling pa- 
rameters are related to each other by To — V{d'^/a?); 
high coupling [Te, 3> 1) for the classical 2DDS hquid 
emulating the closely spaced EHB is now dictated by the 
ordering F^^ <^ d"^ /a? < 1, a condition that is easily met 
for fixed layer density 71 at sufficiently low temperatures 
in the classical regime and/or d sufficiently small. In the 
present work we focus on the classical 2DDS in the strong 
coupling regime that includes both the dipole liquid and 
solid phases. Our preliminary MD study indicates that 
the classical 2DDS liquid freezes at F^i ^ 70. 

In the low temperature regime the 2DDS becomes a 2D 
bosonic dipole liquid, with superfluid properties. The col- 
lective mode spectrum of this model has been considered 
by a number of investigators (20l - [2^ . Quantum Monte 
Carlo (QMC) simulations carried out by Astrakharchik 
and co-workers [13, Hlj paved the way by generating es- 
sential information about the 2DDS ground-state energy, 
static structure function S{k), and the 2DDS liquid-solid 
phase transition, which they estimate to occur at ro ^ 
30. Invoking the zero temperature Feynman Ansatz 
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they established an upper-bound estimate of the col- 
lective mode dispersion with the input of their MC- 
generated S{k) data. Subsequently, Mazzanti and co- 
workers [11] generated collective mode spectra based on 
the more sophisticated correlated basis functions (CBF) 
formalism [26l [27j. using S{k) as an input. 



In contrast to the above studies that require addi- 
tional assumptions to make it possible to infer infor- 
mation about the excitation spectrum from static struc- 
ture function data, an entirely different approach, namely 
classical MD simulations can provide direct insight into 
the structure of the collective modes. We argue that such 
information, even though based on classical dynamics, is 
pertinent both to the classical and to the quantum do- 
mains: it is expected that the nature of the collective 
excitations would not be all that different in the two do- 
mains. This is borne out by our preliminary MD studies 
[13, [13 of the longitudinal collective mode dispersion in 
the strongly coupled classical 2DDS liquid: we have ob- 
served that the classical dispersion reproduces the qual- 
itative features of the collective mode dispersion calcu- 
lated from the Feynman formula ([T]) ^,[ll|, Fig.[T]makes 
this point clear. More precisely, the dispersion curve of 
the classical 2DDS liquid, as generated by our MD simu- 
lations, falls in the narrow band of dispersion curves for 
the 2D bosonic dipole system: the band is bounded from 
above by the Feynman Ansatz (Eq. ([T])) and from below 
by the CBF calculated dispersion. 
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FIG. 1. Longitudinal collective mode dispersion curves gen- 
erated from classical molecular dynamics (MD) simulations 
(squares) at Fd ~ 28 and from the Ref. [iO, [13| zero- 
temperature quantum Monte Carlo (QMC) simulations (dots) 
at td = 28.4. cju = \J2-Kn\i?' jTna? is a characteristic dipole 
oscillation frequency. 

The remarkable likeness of the strongly coupled 
bosonic 2DDS dispersion to the excitation spectrum in 
the superfluid phase of "'He [1^ has been noted by 
Kalnian et al. [24]. The most apparent similarity is in 
the formation of the "roton minimum" , preceded by the 
"niaxon" maximum (see Fig. [1}. Contrasting, however, 
with liquid ^He, the unique feature of the 2DDS resides in 
the simplicity of the purely repulsive inter-particle inter- 
action and the possibility of tuning the coupling strength 
by changing system parameters such as density and tem- 
perature. This makes it possible to analyze the longi- 
tudinal collective mode dispersion and its evolution as a 
function of the F^i and ro coupling parameters defined 
above and to relate them to the experimentally observed 
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phonon/maxon/roton dispersions in ^He. MD simulation 
studies to this effect have been recently carried out by us 
to demonstrate that the emergence of the roton mini- 
mum in bosonic superfluids is, in fact, a consequence of 
strong particle correlations and, as such, is basically a 
classical effect [24| . This finding gives further impetus to 
the present more detailed analysis of the collective mode 
structure of the 2DDS. 

It is in the light of these observations that we have un- 
dertaken a detailed combined analytic/computer simula- 
tion study of the collective mode behavior of the classical 
2DDS liquid and solid phases, whose results are presented 
in this paper. 

The theoretical calculation of the collective mode dis- 
persion in the liquid phase is to be carried out us- 
ing the well-tested quasi-localized charge approximation 
(QLCA) [29, 30] with the input of MD-generated pair dis- 
tribution function g(r) data. The observation that serves 
as the basis for the QLCA is that the dominating feature 
of the physical state of a classical dipolar liquid with 
coupling parameter F^i >> 1 is the quasi-localization of 
the point dipoles. The ensuing model closely resembles 
a disordered solid where the dipoles occupy randomly 
located sites and undergo small-amplitude oscillations 
about them. However, the site positions also change 
and a continuous rearrangement of the underlying quasi- 
equilibrium configuration takes place. Inherent in the 
model theory is the assumption that the two time scales 
are well separated, and that it is sufficient to consider 
the time average (converted into ensemble average) of the 
drifting quasi-equilibrium configuration. The calculation 
of the lattice phonons in the solid phase is to be carried 
out using the well-known harmonic approximation. 

MD simulations have been pursued along the lines used 
for simulating other strongly interacting classical systems 
[31I [3^ . The topical review article fSr] provides an in- 
depth description of the computational methodology fol- 
lowed in the present work. 

In the present work, we study the collective excita- 
tions in the 2DDS liquid and solid over a wide range 
of wave numbers. In the liquid phase, both the longi- 
tudinal and transverse modes have been identified. In 
the solid phase, the two modes persist, but the system 
is anisotropic and the polarizations become wave vector- 
and angle-dependent. In our earlier work (ssj . we have 
already analyzed the fca — ?> longitudinal acoustic be- 
havior. Here we go further: Our analysis extends beyond 
the ka > 2 domain where a well-defined maximum of 
the dispersion curve develops followed by a deep mini- 
mum. These extrema are analogous to what has become 
known as the "niaxon" and "roton" excitations in the 
mode spectrum of liquid helium. For this reason, we will 
use the same terminology throughout this paper. Our 
recent MD study of collective modes in the 2DDS solid 
2J| indicate the existence of remarkable roton-|-maxon, 
maxon-|-maxon, roton-|-roton combination frequencies at 
wave numbers in the vicinities of the maxon and the ro- 
ton, again, in analogy with what has been theoretically 



predicted [2a, [SJ-jS^ in relation to liquid ''^He spectra. In 
the present work, we continue this line of investigation 
for the solid and for the liquid phase as well. In rela- 
tion to the phonon spectra we are able to make contact 
with the elastic theory of solids and identify the principal 
elastic constants. 

The organization of the paper is as follows: In Sec. 
II, we analyze the collective mode dispersion. We gen- 
erate longitudinal {L) and transverse (T) spectra from 
our MD simulations of the 2DDS liquid and solid phases. 
The ensuing MD dispersion curves, constructed from the 
peaks of the spectra, serve as standards for comparison 
with the theoretical QLCA L and T oscillation frequen- 
cies displayed in the same section. In Sec. HI, we con- 
sider the solid phase and apply the harmonic approxi- 
mation (HA) to the calculation of the longitudinal and 
transverse acoustic phase velocities, (from which we infer 
the elastic constants of the crystal), and the full phonon 
dispersions for the 2DDS. Conclusions are drawn in Sec. 
IV. 



II. CLASSICAL DIPOLE LIQUID: COMBINED 
MD/QLCA ANALYSIS 

We turn now to the analysis of the longitudinal and 
transverse collective modes in the classical 2DDS liquid. 
This analysis has been carried out using the QLCA the- 
ory combined with MD simulations. 

The successful application of the QLCA to the calcula- 
tion of collective mode dispersion in a variety of strongly 
coupled charged particle systems has been well docu- 
mented over the past two decades. Its recent application 
to the 2DDS liquid has resulted in an accurate descrip- 
tion of the longitudinal collective mode dispersion in the 
acoustic domain, as borne out by tabulated comparisons 
with MD and thermodynamic sound speed data [sS]. 

Our program of MD simulation of the dynamics of 
the classical 2DDS liquid has been carried out along the 
lines used for simulating other strongly interacting clas- 
sical systems [3l|. The present MD simulation involves 
4200 particles; information about the collective modes 
and their dispersion is obtained from the Fourier trans- 
form of the correlation spectra of the microscopic density 



nk 



(2) 



yielding the dynamical structure function 

5(fc, Lj) = {1/27TN) lim (l/AT) \nk{Lu)\^ , (3) 

A— >(X) 

where AT is the duration of the data recording period. 
Similarly, the spectra of the longitudinal and transverse 
current fluctuations, L{k,uj) and T{k,uj), can be ob- 
tained from Fourier analysis of the microscopic currents 



^fe(i) = y^,V]x{t) exp[ikxj{t)] 



(4) 
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and 



'''hit) = ^ Vjy{t) exp[ikxj{t)] 



(5) 



we assume that k is directed along the x axis (the sys- 
tem is isotropic). The collective modes are identified as 
peaks in the fluctuation spectra. The widths of the peaks 
provide additional information about the lifetimes of the 
excitations: narrow peaks correspond to longer lifetimes, 
while broad features indicate short-lived excitations. 

Before proceeding with the analysis, a recap of our 
recent findings psl Issj for the fca — )■ regime is in order. 

The crucial observation concerning the collective be- 
havior of the 2DDS is that average the Hartree field does 
not exist, since 



{(f>Dir))H =n (fr4)D{r) 



(6) 



is unbounded; therefore, the Fourier transform of the 
dipole potential = /i^/r^ does not exist, implying 

that the 2DDS can have no RPA limit. Comparing this 
situation with the case of a 2D Coulomb system where 
either the routine observation 



i:j^(fc) oc fc^(/)couL(fc), (/)couL(fc) oc 1/fc, 
or a simple scaling argument for the integral 



0(k) — J d^r 4>{r) exp(zk • r) 



provide the correct 



uj(k) cx Vk 



(7) 



(8) 



(9) 



RPA plasmon behavior, a similar line of reasoning in the 
present case would lead, via the (faulty) relation 



ipoik) oc -k^ (t>GOVi.{k) cx fc. 



(10) 



to the incorrect RPA longitudinal collective mode fre- 
quency as [III 



(11) 



The correct approach, as established in Ref. [33[, leads 
to a correlation controlled long- wavelength acoustic dis- 
persion 



i[k) cx k. 



(12) 



of the strongly coupled 2DDS liquid both in the classical 
and in the zero-temperature quantum domains 28, 33]. 
This behavior is also similar to that of the EHB liquid 
[39| , with phase velocities virtually identical in the two 
cases. The precise numerical value of the acoustic veloc- 
ity is governed by the average potential 



(r)) = n / d^r (t)D{r)g{r), 



(13) 



rather than by the Hartree potential; g{r) is the 
statistics- and coupling-dependent equilibrium pair dis- 
tribution function. We found that the theoretical QLCA 
values of the longitudinal acoustic phase velocity, cal- 
culated from the Eq. ([13]) integral agree reasonably well 
with the MD simulation values over a wide range of classi- 
cal liquid phase coupling strengths; the disparity between 
the theoretical acoustic phase velocities and thermody- 
namic sound speeds, though somewhat larger, is still only 
about 5.3%. The longitudinal acoustic velocity decreases 
with increasing coupling parameter and assumes its low- 
est value in the solid phase. 

The details of the derivation of the longitudinal and 
transverse dispersions for arbitrary k values based on 
QLCA theory culminating in the 2DDS dispersion rela- 
tion below are given in [33j]. Here it suffices to take 
Eq. (IT4t , together with the QLCA dynamical tensor (|15l) , 
as the appropriate starting point for the present study: 



a; 



'6, 



C^.Ml =0, 



(14) 



C^,(k) =-- [ A\ g(r)[exp(zk. r) - l\d^,dMr) (15) 
m J 

3n^2 



d r — g(r)[exp(ik • r) — 1] 



The longitudinal (L) and transverse (T) oscillation fre- 
quencies are readily calculated to be 



wi(fc) = CUk) (16) 

= -u:l / df -g{f) [3 - iUkf) + ^Mkf)]] 
Jo f 

u:l{k) = Crik) (17) 

3 f°° 1 
= -u:l / df -g{f) [3 - 3Jo(fcf) - 5J2(fcf)]] , 



— 27rn/i^/ma'^, f = r/a, and k = ka. Note that the 
usual interpretation of Eqs. (|T6l) and (|T7| would require 
splitting g{r) as 



g{r) = 1 + h{r) 



(18) 



and identifying the frequency coming from the "1" piece 
in the integral as the RPA, while h{r) as the correlational 
contribution. The discussion presented above shows the 
fallacy of this reasoning for the 2DDS. For Tn fixed, dis- 
persion curves can be generated from (|16p and ()17p with 
the input of g{r) pair distribution function data obtained 
from our MD computer simulations. Sample data for the 
latter are displayed in Fig. [2] 

The wave number regimes of special interest are (i) 
the long-wavelength (fca — > 0) acoustic regime, where 
both longitudinal and transverse acoustic modes develop. 



5 



3 
2.5 
2 

ra 1.5 
1 

0.5 




1 






rD = 7 

- 15 

rD = 30 
Pr, - fin 
























i 

_ I 


\| 




















- 
















1 











6 

r/a 



10 



12 



FIG. 2. MD pair distribution function for To = 7, 15, 30, 
and 60. 



(ii) the finite wave number regime spanning the maxon- 
roton portion of the longitudinal dispersion curve, and 

(iii) the high-Zc domain where the dispersion is dominated 
by single particle excitations. 

Addressing first the long-wavelength (fca — !■ 0) regime, 
Eqs. (dH) and (H?]) simplify to 



0) = |j(rz,)..!,P, 



0) 
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2 0D(a) 



dr^g{r). 



(19) 



(20) 



(21) 



The values of JiTjj) calculated from (|2ip with the in- 
put of MD-generated pair distribution function data and 
the ensuing values of the QLCA acoustic velocities cal- 
culated from Eqs. ([T^ and (PU)) are displayed in Table I 
along with the corresponding MD acoustic phase veloc- 
ities and the thermodynamic sound speeds (s™) calcu- 
lated from the isothermal cornpressibility using the clas- 
sical 2DDS equation of state [3J| . Using the MD data as 



TABLE I. 2D-point dipole liquid: QLCA (sg'^j?^), MD {s^°), 

and thermodynamic (s"'"'^) sound speeds as functions of the 
classical coupling parameter Fd- The entries in columns 2-5 
are quoted from [3J| for To = 10 — 40, too = ^2im^'^ /ma^ . 
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10 


0.8847 


1.351 


1.312 


1.282 


0.4073 




20 


0.8504 


1.324 


1.276 


1.257 


0.3992 




30 


0.8370 


1.314 


1.246 


1.247 


0.3962 


0.38 


40 


0.8295 


1.308 


1.258 


1.242 


0.3944 


0.38 


60 


0.8208 


1.301 


1.256 


1.234 


0.3923 


0.37 



a reference, the discrepancy between longitudinal QLCA 
and MD sound speeds ranges from 2.89% at Fd = 10 to 
3.46% at F^i = 60. The entries for the MD transverse 
acoustic speeds in column 7 are calculated as the differen- 
tial slope s — dio/dk in the linear regime at fca « 1. Such 
data for Fij < 20 are not available because in this weaker 
coupling regime the shear wave is too heavily damped, 
similarly to the situation for the 2D0CP [40]. The dis- 
crepancy between the transverse QLCA and MD acous- 
tic velocities ranges from 4.09% at Fj^ = 30 to 5.68% at 
F^ = 60. 
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FIG. 3. (color online) Longitudinal (L) current fluctuation 
spectra for Fd ~ 60; a = l/^/nn is the 2D Wigner-Seitz 
radius and cjd ~ \J2-Knii?- jmcfi is a characteristic dipole os- 
cillation frequency. The thick arrows point in the direction of 
increasing fco. 



In Table U note that multiplication by -v/STd con- 
verts auJD units into l/y/(3m thermal velocity units; 
u]D = ^j2TTn^j?Jrria? is a characteristic dipole oscilla- 
tion frequency. In these latter units, the thermodynamic 
speed decreases from a value of 13.52 at = 60 to unity 
at Tjj = 0. We observe that the MD and QLCA longi- 
tudinal phase velocities in Table |T] are somewhat higher 
than the thermodynamic speeds. The analysis of the 
2DDS crystal in Sec. 3 may shed some light on the ori- 
gin of this difference. 

Turning next to the finite ka domain, we have gener- 
ated MD spectra for a wide range oiTo values. Repre- 
sentative L(k,u;) and T(k,uj) spectra at Tjj — 60 are 
displayed in Figs. [3] and SI In Fig. [3] there are four 
ka domains between ka — 0.05 and ka = 5.66. Look- 
ing at Fig. [3^ we observe that the spectral peak shifts 
to higher frequencies with increasing ka values, reach- 
ing the maxon frequency lum ~ 1.2a;£) around fca ~ 2 
(Fig. [SJj). For ka values in the sub-interval [1.48,2.64] 
depicted in Fig. |3Jd, we observe two well-defined clus- 
ters of spectral peaks: the lower frequency cluster rep- 
resents the domain of the maxon; within the higher fre- 
quency cluster, the broad peaks indicate the emergence of 
a faint maxon-maxon (M-l-M) harmonic. For ka values in 
the sub-interval [2.58,3.90] depicted in Fig. [31:, the lower 
frequency peaks shift to the domain of the roton mini- 
mum which assumes its lowest frequency value at around 
ka w 3.7; we note the persistence of the maxon-|-maxon 
harmonic over this entire ka sub-interval. At the higher 
ka values depicted in Fig. [5Ji, the cluster of lower fre- 
quency spectral peaks indicates the existence of a second 
maximum at around ka 5.4. 




FIG. 4. (color online) Transverse (T) current fluctuation spec- 
tra for Td = 60. 



The resulting MD longitudinal and transverse disper- 
sion curves obtained from the peaks of the L and T spec- 
tra, respectively, are displayed in Figs. [5] and [6l 

Addressing first the longitudinal mode dispersion in 
Fig. [5l we observe the progressive deepening of the ro- 
ton minimum with increasing Td. The broadening of 
the spectral peaks becomes more and more pronounced 
with decreasing coupling strength; that is, the lifetime 
of the collective mode decreases with increasing temper- 
ature. The simulations indicate that, for wave numbers 
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FIG. 5. Longitudinal dispersion curves for Vd ~ 7,15,30, 
and 60. The vertical lines indicate the widths of the spectra. 
The shaded region in (d) should be regarded as a ghost-like 
dispersion depicting the emergence of a faint maxon-maxon 
(M-l-M) harmonic. 



< fca < 2.5, spanning the entire acoustic domain and 
extending somewhat beyond the maxon peak, the longi- 
tudinal mode is fairly robust. At higher wave numbers, 
the mode is strongly damped, but can still be viable over 
the approximate domain 2.5 < fca < 5 so long as the 
2DDS liquid is in the Fd > 30 strong coupling regime. 
At F/3 = 60, we observe from Fig. [5] that the location of 
the roton minimum in the classical 2DDS liquid is close to 
the location of its counterpart in the 2D bosonic dipole 
liquid at zero temperature [13, Hlj. In the Conlusions 
we further discuss the relevance of the observed roton- 
maxon behavior to the dispersion characteristics of the 
2D bosonic dipole liquid. 

In the high-fca domain, we observe that the longitu- 
dinal dispersion eventually becomes dominated by the 
Bohm-Gross (BG) oscillation frequency 

, 3 3 kcL , . 

characteristic of single-particle excitations in the classical 
2DDS, and in strict compliance with the third- frequency- 
moment sum rule for dipole systems 41| . The switchover 
ka value for the single particle beavior occurs between 
ka = b (for Td = 15) and ka = 8.5 (for F^, = 60). 
Again, making contact with the 2D bosonic dipole liq- 
uid, it is of some interest to compare this a;(fc — oo) oc k 
asymptotic behavior with the one predicted by the Bo- 
golyubov theory for the quasicondensate 2DDS [H, [3l| . 
In this latter the high-fca behavior sets on around simi- 
lar fc- values ^] and is also dominated by single-particle 
excitations, which, however, originate, in sharp contrast 
to the classical system, from the zero momentum ground 
state. Hence the difference in the asymptote, which in 
the zero temperature limit is aj(fc — >■ c») oc hk'^ /2m. 
The shaded region in Fig.lSfd) is seen as the emergence 
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of the faint maxon-maxon (M+M) harmonic detected in 
the Fig.[3]spectra. More will be said about this harmonic, 
along with other combination frequencies that emerge 
only in the 2DDS lattice, in Sec. III. 




FIG. 6. Transverse shear mode dispersion curve for Td ~ 60. 
Note the wave number cutoff fc*a as a function of Fd marking 
the ka value where uj{ka) — 0. 

We next address the transverse shear mode dispersion 
displayed in Fig. [6 Similarly to the 2D one-coniponent 
plasma (2D0CP) i^l and to complex plasmas (l^ - li^ . 
shear waves in the strongly coupled 2DDS liquid are 
strongly damped, even at F^i = 60, as evidenced by the 
sizeable line widths shown in Fig.jB] for Fd < 30, the MD 
simulations indicate that the shear waves are too heavily 
damped to be viable. In the coupling regime where they 
are viable. F ig, p and Fig.|9]below show that, similarly to 
the 2D0CP [loTand to 2D Yukawa plasmas (2DYP) [il- 
2DDS liquid-phase shear waves cease to exist below 
a critical finite F^i-dependent wave number cutoff, fc*a, 
marking the ka value where u){ka) = 0. As expected, 
fc*a ^- as Fz5 T^owd ^^^^ below. 

Turning now to the QLCA description of the mode dis- 
persion at finite wavenumbcrs, the straightforward calcu- 
lation of the QLCA oscillation frequencies ([T6l) and ([T7| , 
with the input of MD-generated g{r) data, results in the 
longitudinal and transverse dispersion curves displayed 
in Fig.mfor F^, = 7, 15, 30, and 60. Figs. Hand [9] pro- 
vide a comparison between the QLCA dispersion curves 
and the MD data. 

Figure [5] shows very good quantitative agreement be- 
tween theory and simulation up to /ca ~ 2.5. As ex- 
pected, the ka range of agreement increases with in- 
creasing coupling, since the QLCA theory is, after all, 
premised to be a strong coupling theory. For increas- 
ing fca-values, we can observe the evolution of the roton 
minimum with increasing F/j. According to the QLCA 
description, the position of the roton minimum shifts to 
lower ka values and deepens more and more with increas- 
ing coupling strength. The MD data, however, indicate 
that the ka ~ 3.7 position of the roton minimum re- 
mains more or less the same as it progressively deepens 
with increasing coupling strength. While the QLCA does 
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FIG. 8. MD and QLCA longitudinal dispersion curves for 
Fd = 7, 15, 30, and 60. 



capture the qualitative features of the roton portion of 
the 2DDS dispersion curve for F^i > 30, quantitative 
agreement between theory and simulation for fca > 3 is 
less satisfactory, even in this high coupling regime of the 
liquid phase. 

Turning now to the transverse shear mode. Fig. [3] 
shows only fair agreement between the QLCA and MD 
dispersion curves. Similarly to what has been reported 
for 2D Yukawa liquid plasmas [i^, the MD data, an- 
chored by the (fc*a, a; = 0) point, lie to the right of the 
QLCA dispersion curve and run parallel to it in the linear 
regime. This discrepancy is due to the QLCA's inability 
to account for diffusional and other damping effects that 
preclude the existence of long- wavelength shear waves in 
the liquid phase In the case of 2D Yukawa liquids 
(isf , the discrepancy was removed by introducing a phe- 
nomenological diffusional damping time into the QLCA 
formalism. We expect that a similar such modification in 
the 2DDS QLCA formalism would bring about the same 
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FIG. 9. MD and QLCA transverse dispersion curves for Yd 
60. 



outcome. 



III. DIPOLE LATTICE: HARMONIC 
APPROXIMATION AND MD SIMULATIONS 

According to our MD simulation the 2DDS liquid 
freezes at about ~ 70. In this Section we turn to the 
analysis of the collective modes (phonon dispersion) of 
the lattice phase. As in all systems with an isotropic cen- 
tral force interaction, the 2DDS crystallizes in a triangu- 
lar (hexagonal) lattice. We calculate the phonon disper- 
sion in the harmonic approximation (HA). The starting 
point is the dispersion relation ([1]) with the liquid-phase 
dj^namical tensor (fT5|) replaced by its HA counterpart 



C^,(k) = -1 ^[exp(ik • rO - l]d^,dAD{n) (23) 

i 

= ^^[exp(zk.r,)-l]4 

TO ^ — ' 1^ 



where k • r.^ = kriC0s{9i — (p);9i is the angle of the 
vector and (f is the propagation angle measured from the 
axis pointing toward the nearest neighbor. 

At long wavelengths, the hexagonal lattice is isotropic 
so that the direction of k is arbitrary. So for conve- 
nience, one can choose k to be along one of the two crys- 
tallographic axes, say, (/9=0". In this case Eq. (|23)) then 
simplifies to 



longitudinal and transverse oscillation frequencies then 
follow from the dispersion relation ([Ti)) and ([M)) : 



a;|(fc 0) = CL(fc ^ 0) (25) 
_^t^j2cosH,^ [l-5cos2 0,] , 



CJ^(fc 



0) 



CT(fc^O) 



(26) 



cos^ 0i— [l - 5sin^ 9i] . 



The hexagonal lattice can be decomposed into a se- 
quence of concentric hexagons tilted away from the crys- 
tallographic axes and of increasing size. Thus the lattice 
sum (1251) can be decomposed into an ri-dependent part 
and an angular part; this latter can be summed over the 
vertices of each of the hexagons. This summations is fa- 
cilitated by observing that for any arbitrary angle a, the 
sum J2n=o cos^(a -I- n7r/3) = 1/2 is independent of a. It 
therefore suffices to sum the 9i portion of the summands 
in (1251) and (j26p around a single hexagon and divide by 
six. This results in 



iolik ^ 0) - CL(fc ^ 0) = —Mujla^k 



ojl{k ^ 0) = Crik 



0) = -MioWk', 



(27) 



(28) 



where M = ^^"^ lattice sum over the trian- 

gular lattice with = r^/a. The longitudinal oscilla- 
tion frequency (j27p has been reported in [33] and is dis- 
played here alongside the new transverse result (l28l) for 
comparison. In effect, the lattice sum M/2 replaces its 
liquid-phase counterpart integral J{Tu) in Eqs. ([T51) and 
(PU]) . The value of M has been calculated by a number of 
workers [46h49|] with slightly different results; the most 
recent semi-analytical calculation due to Rozenbaum [i^ 
is quoted here as M = 1.642. Our own lattice sum com- 
putation for the 2DDS crystal involving 1.9 x 10^ particles 
provides M = 1.597 From Eqs. ^ and (HH), the 
corresponding L and T sound speeds are 



(29) 



C^.(fc^O) 



21,2 



cos^ 9i 



SuL' - 5 



(24) 

Moreover, for the isotropic system, the dynamical tensor 
(I24p is diagonal, implying that the eigenvectors are par- 
allel and perpendicular to k, as in the liquid. The purely 



St = 0.387wDa. 



(30) 



As expected, these values are slightly lower than their 
respective = 60 liquid phase counterpart entries in 
Table m Note the agreement between (P^)) and the mea- 
sured value 1.2836ati;£i- The thermodynamic sound speed 
s"""^, though not physically meaningful for a lattice, is 
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still of interest from the point of view of providing an 
estimate of the liquid phase thermodynamic sound speed 
in the large-Fu limit. Table U indicates that the ther- 
modynamic sound speed decreases monotonically with 
increasing F^i. One therefore expects that s™ assumes 
its minimum value when Td oo. 

For an isotropic elastic medium the long wavelength 
{k — ^ 0) behavior of the phonons can be described in 
terms of longitudinal and transverse elastic waves, rather 
than in terms of the thermodynamic sound speed, s™, 
with phase velocities expressed in terms of the elastic 
constants K and G. K — —V{dP/dV) is the bulk mod- 
ulus, G is the shear modulus. In 2D the corresponding 
velocities are 




Ik + g 




(31) 



with p = mn being the mass density of the solid. It 
is important to note the marked difference between the 
thermodynamic sound veloicity s"^^ and the elasic phase 
velocities sl and st as enunciated by (|31|) . There is a 
small, but discernible, difference between the s™ and sl 
and St, 

.TH ^ y^TI^ (32) 

readily follows from (1511) . From and pOI) . the ratio 
s^l s\ — 1/11 then provides 

s™ = sl/1.049 = 1.223aa;£, (33) 

for the 2DDS hexagonal lattice. This value is in keeping 
with the trend in Table ID 

We address next the lattice phonon dispersion at finite 
wavenumbers obtained by solving the dispersion relation 
(fn)) with the input of the dynamical tensor . The re- 
sults very closely resemble the dispersion of phonon spec- 
trum of another known 2D swtem, the 2D Yukawa tri- 
angular (hexagonal) lattice [3l|. In particular, the value 
1/11 of the sound speed ratio s\/s\ turns out to be the 
same as the sound speed ratio for the 2D Yukawa crys- 
tal [131 with a screening parameter na = 1.05. Phonon 
dispersion curves displayed in Fig. [10] for four propaga- 
tion angles: ip = 0°, 10", 20^ 30° {ip = O" and ip = 30° 
are the crystallographic axes). The angle 0, indicated in 
the right panels of Fig. IH is the polarization angle mea- 
sured with respect to the propagation vector k; the mode 
polarizations are purely longitudinal (Q — 0°) or purely 
transverse (0 = 90°) for propagation along the = 0° 
and ip = 30° crystallographic axes only. Otherwise, the 
polarizations are mixed as shown in the ip = 10° and 
If = 20° right panels of Fig. [S] The dispersion curves 



are periodic in fc; the period is the reciprocal lattice con- 
stant only for propagation along a crystallographic axis; 
for intermediate angles, the much longer period is given 
by formulas (45) and (46) in [3T'|. 
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FIG. 10. Electric dipole hexagonal lattice eigenmodes and 
corresponding mode polarizations for angular directions (with 
respect to the nearest neighbor direction) <^ = — 30° in lO'' 
steps. The measured sound velocity sl = 1.284aa;D. 

For a macroscopically disordered lattice one can view 
the system as an aggregate of locally ordered domains 
with randomly distributed crystallographic axes. One 
may also contemplate this as an alternate model for the 
strongly coupled liquid. We can then seek the simi- 
larity to the liquid-phase dispersion through a suitably 
angle-averaged dispersion of the lattice. The procedure, 
which parallels that of Ref. [3l|, consists in projecting out 
the longitudinal and transverse components of the lat- 
tice eigenmodes and comparing their respective angular 
averages with the longitudinal and transverse modes in 
the liquid phase. The calculation of the angle-averaged 
longitudinal and transverse lattice dispersions proceeds 
according to the following prescription: First, we com- 
pute the lattice eigenmodes for k vectors with angular 
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directions between ip = and (p = 30° in one degree 
steps. We next perform the longitudinal and transverse 
projections based on the normal mode data with 



, ,2 _ 



(34) 



k • Gi 



k • §2 



^2' 



where k, ei and §2 are unit vectors parallel to the wave 
vector (k) and normal-mode eigenvectors (ei and 62), 
resp ectively; uji and 0^2 are the normal-mode frequencies 
|51| . We can then simply average the frequency values 
belonging to equal k values for the longitudinal and trans- 
verse dispersions separately. 
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FIG. 11. 2DDS lattice dispersion. Angle- averaged longitudi- 
nal and transverse lattice modes shown with liquid phase MD 
data and the QLCA dispersion curves at Vo = 60. 

The results are displayed in Fig. [TTJ Addressing first 
the agreement between the two approaches and the MD 
we see that the QLCA provides a superior descrition for 
< /ca < 2.8 and somewhat less satisfactory quantita- 
tive agreement thereafter. In contrast the lattice average 
dispersion, while it shows excellent agreement with the 
MD data over the narrower interval < /ca < 1, com- 
pares somewhat less favorably with the MD data than 
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FIG. 12. (color online) Illustration of the appearance of com- 
bination frequencies in the longitudinal current-current cor- 
relation functions L{k = kx,ijj), measured along the direction 
of the nearest neighbor in a hexagonal lattice simulated at 
F_D = 1000. The center color map of the L{k = kx,oj) fluc- 
tuation spectra shows the strong primary dispersion and the 
ghosts of the combination frequencies, thin white lines are the 
dispersion curves from both L{k = kx,u)) and T{k — kx,u)). 
The four small panels show vertical cross-sections taken for 
the selected group of wave numbers, where the appearance of 
the combination frequencies is the most manifest. 



the QLCA thereafter, but becomes superior to the QLCA 
in the neighborhood of the high-fc roton region. As to the 
transverse collective mode dispersion shown in Fig. [TTb . 
we observe the same trend. The discrepancy between the 
QLCA dispersion and the MD data was discussed earlier 
in connection with the finite ka cutoff. 

In Figs. [12] and [T21 we have displayed a series of 
L[k = kx^Lo) and L(k — ky^Lo) spectra generated over 
a wide range of ka values dXTn — 1000. Our findings 
reaffirm the observations from our previous simulations 
[23 | carried out at Vu = 500 and over a narrower range of 
wave numbers: The accumulation of weight in the vicin- 
ity of three combination frequencies, labeled as "1" , "2" 
and "3" is clearly visible. The emergence of harmonics 
appears to be a general feature of Yukawa, Coulomb, and 
other types of similar interaction (52j as well, although 
the quantitative relationships between the amplitudes of 
the different combination frequencies should be sensitive 
to the actual interaction potential. 
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FIG. 13. (color online) Same as Fig. [12] for wave-numbers 
perpendicular to the nearest neighbor direction, L{k = ky,u)). 



The identification of tlie "1" , "2" , "3" peaks witfi tlie 
roton-roton (R+R) , maxon-roton (M+R and perhaps M- 
R), maxon-maxon (M+M) combinations, respectively, is 
tempting, but the full understanding of this process re- 
quires more study. It should be realized that due to the 
anisotropy of the medium, the formation of combination 
frequencies in a lattice is a much more intricate matter 
than in an isotropic liquid. First, since the positions of 
the minima and maxima of the dispersion curves are an- 
gle dependent, the notions of roton frequency and maxon 
frequency become ill-defined and the strengths and po- 
sitions of the combination frequencies become also an- 
gle dependent; this is well illustrated by the comparison 
of Figs. [12] and [T3] where modes propagating along the 
two different crystallographic axes are compared. Sec- 
ond, most of the weight should now originate from the 
vicinity of the local extrema of the two-dimensional dis- 
persion surface, rather than from the extrema of the se- 
lected dispersion curves; but the nature and locations 
of these former have not yet been explored. Third, the 
appropriate matching of the fc-vectors may become a deli- 
cate issue, and one would expect that directions in which 
the highest amplitude harmonics appear may be quite 
different from the directions associated with their parent 
frequencies. In fact, as pointed out by [13], it is the very 
close proximity of the fc = region where the largest 
possible number of combination k-vectors could end up. 



thus leading to the expectation that high- amplitude com- 
bination frequencies would appear in this domain. It is 
therefore remarkable that our data do not indicate the 
presence of any weight in the fc = region: whether this 
is due to the limitations of the simulation protocol or a 
real effect, still remains to be understood. 



IV. CONCLUSIONS 

In this paper we have carried out a combined an- 
alytic/molecular dynamics (MD) study of the dynam- 
ics and collective mode dispersion of a classical two- 
dimensional dipole system (2DDS) in its strongly cou- 
pled liquid and solid phases. The analytical methodol- 
ogy is based on the quasi-localized charge approxima- 
tion (QLCA) for the description of the liquid-phase dis- 
persion, and on the companion harmonic approximation 
(HA) for the description of the lattice phonons. As we 
stated in the Introduction to this paper, we had two ob- 
jectives in mind: The first objective was to understand 
how the dynamics of classical many-body system with 
a 1/r'^ interaction potential are related to the dynam- 
ics of similar systems, such as the 2D OCP and the 2D 
Yukawa plasma. The second objective was to relate the 
collective mode behavior of the classical 2DDS to that of 
its quantum counterpart, the 2D bosonic dipole system. 

As to the first objective, in contrast to the 2D OCP 
which features an cj(fc — > 0) oc ^/k longitudinal mode, the 
2DDS exhibits an u;(/c — >■ 0) oc fc acoustic excitation. In 
this sense, the behavior of the 2DDS is similar to that of 
the 2D Yukawa system. The 2DDS, however, is unique in 
that the longitudinal acoustic mode is wholly maintained 
by strong short-range correlations. 

Similarly to the 2D0CP [2^ and complex plasmas 
our MD simulations reveal that transverse shear 
waves propagating in the strongly coupled 2DDS liquid 
are strongly damped, especially for coupling strengths 
To < 60. Moreover, the waves cease to exist below a 
critical wavenumber cutoff, marking the point where the 
excitation frequency becomes identically zero. As ex- 
pected, this critical wave number decreases with increas- 
ing coupling strength, tending to zero when the 2DDS 
liquid freezes. 

Using the harmonic approximation, we have calculated 
the longitudinal and transverse phonon dispersions of the 
crystalline solid phase. The morphology of the lattice 
phonons closely emulates that of the 2D Yukawa lattice; 
this is not surprising since the morphology appears to be 
dictated primarily by a common hexagonal lattice struc- 
ture rather than by the nature of the interaction. For 
the macroscopically disordered lattice, where the system 
can be viewed as an aggregate of locally ordered domains 
whose symmetry axes are randomly distributed, we have 
calculated the angle-averaged lattice dispersion curves. 
One may also contemplate this approach as a model al- 
ternate to the QLCA for the strongly coupled liquid. The 
comparison of the two of them with the MD results show 
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reasonable agreement up to high ka values: the QLCA 
seems to perform better in the lower, while the averaged 
lattice model in the higher ka regions. 

Addressing the second objective, we note that the main 
interest in this regard lies in the understanding of the 
formation of the so-called maxon-roton structure of the 
dispersion in the 2D bosonic dipole liquid, which is also 
shared by liquid '^He. The terms "maxon" and "roton" 
serve as convenient labels for the first maximum followed 
by a deep minimum of a longitudinal dispersion curve. 
We have found that in this domain and in the strong 
coupling regime, there is a close affinity between the dis- 
persion curves of the classical 2DDS liquid and the low- 
temperature bosonic dipole liquid. More precisely, in a 
recent work [2^ we have shown that the dispersion curve 
of the classical 2DDS liquid, as generated by our MD 
simulations, falls in the narrow band of dispersion curves 
for the 2D bosonic dipole system: the band is bounded 
from above by the Feynman Ansatz (Eq. ([T])) and from 
below by the CBF calculated dispersion. 

Thus, we believe that a detailed analysis of the col- 
lective mode behavior of the classical 2DDS liquid will 
reflect - except for the damping - the collective mode 
behavior of its quantum counterpart. As to the damp- 
ing, the high-Zca modes, in general, are strongly damped 
in a classical liquid to the extent that they become unob- 
servable. This seems to be the case in the 2D OCP and 
Yukawa plasmas. Here, we have found that the expected 
strong damping prevails in the 2DDS liquid as well. By 
contrast, the collective modes in the 2D quantum system 
in its superfluid phase are virtually undamped. Thus, we 
expect that the collective mode features that we analyze 
in the classical domain will correspond to the observ- 
able features in the collective modes of the 2D bosonic 
dipole system in its superfluid phase. We have considered 
three wave number domains for the analysis of the collec- 
tive modes in the liquid phase: (i) the long-wavelength 
{ka — > 0) domain, where both longitudinal and trans- 
verse acoustic modes develop, (ii) the finite wave number 
domain spanning the maxon-roton portion of the longitu- 
dinal dispersion curve, and (iii) the high- A; domain where 
the dispersion is dominated by single-particle excitations. 
Our findings from the MD simulations are as follows: For 
wavenumbers < fca < 2.5 extending somewhat beyond 
the maxon portion of the dispersion curve, the longitu- 
dinal collective excitations are robust over a wide range 
of coupling strengths 7 < Fd < 60. During the pro- 
gressive deepening of the roton minimum with increas- 
ing coupling, its ka ~ 3.7 position remains more or less 
the same. In keeping what we have already stated, this 



wave number value in the classical 2DDS liquid is quite 
close to the location of the roton minimum in the quan- 
tum 2DDS at zero temperature [l3|- For wave numbers 
well above the roton minimum, the dispersion assumes 
the character of single particle excitations. This feature 
is also qualitatively the same - albeit with different k- 
dependence -, as the one predicted and observed in the 
high-fc domain for the 2D cbosonic dipole system. 

Additionally, we observe the existence of a faint 
maxon-maxon (M-l-M) harmonic frequency that persists 
over the broad range of wave numbers 0.05 < A:a < 5. 
This feature can be related to the observation of roton- 
roton, roton-maxon, and maxon-maxon combination fre- 
quencies in the superfluid phase of ''He. Whether the 
absence of the first two combination frequencies is a fea- 
ture that distinguishes the 2DDS from liquid ''He, due 
to the difference in interaction potentials, or these two 
combination frequencies are simply masked by classical 
noise is not clear. More combination frequencies have 
been observed in the crystalline solid phase. However, 
due to the anisotropy of the medium, the formation of 
combination frequencies in a lattice is a much more in- 
tricate matter than in an isotropic liquid. The definite 
identification of the frequencies with the maxon-maxon 
(M-l-M), maxon-roton (M-l-R and perhaps R-M), roton- 
roton (R+R) combinations would still require more study 
and further understanding of this process. 

We have compared the MD data with the theoretical 
analysis based on the QLCA approach. As expected, the 
range of ka values marking agreement between QLCA 
theory and MD data increases with increasing coupling. 
At 0.05 < ka < 2.5 there is near-perfect agreement be- 
tween theory and simulation. While, however, the QLCA 
dispersion closely emulates the qualitative features of the 
dispersion even at higher ka values, in the neighborhood 
of the roton minimum it underestimates the roton fre- 
quency and overestimates its wave number. An improve- 
ment of the QLCA, with the capability to account for the 
strong angular correlations in the strong coupling domain 
may be needed to deliver even better agreement with MD 
observations. 
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